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^ ■ Abstract 

O 
o 

In a previous paper the statistical properties of the hnearized Kolmogorov 
^ ' flow have been studied, using the formalism of fluctuating hydrodynamics. In 

this paper the nonlinear regime is considered, with emphasis on the statistical 
properties of the flow near the first instability. The normal form amplitude 
equation is derived for the case of an incompressible fluid and the velocity 
field is constructed explicitly above (but closed to) the instability. The rel- 
ative simplicity of this flow allows one to analyze the compressible case as 
well. Using a perturbative technique, it is shown that close to the instability 
threshold the stochastic dynamics of the system is governed by two coupled 
non linear Langevin equations in the Fourier space. The solution of these 
equations can be cast into the exponential of a Landau Ginzburg functional, 
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which proves to be identical to the one obtained for the case of the incompress- 
ible fluid. The theoretical predictions are confirmed by numerical simulations 
of the nonlinear fluctuating hydrodynamic equations. 
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I. INTRODUCTION 



A central issue in nonequilibrium statistical physics is the role of fluctuations in the onset 
of hydrodynamical instabilities. From a theoretical point of view one generally relies on 
the Landau-Lifshitz fluctuating hydrodynamics 0, mainly because of its relative simplicity 
as compared to more fundamental approaches ||^^. Fluctuating hydrodynamics has been 
used by various authors to study the statistical properties of simple fluids subjected to 
nonequilibrium constraints, such as temperature gradient or shear (for a review, 

see ref. [|]). Light scattering results, obtained for systems under temperature gradient, have 
shown quantitative agreement with theoretical predictions Quantitative agreement has 
also been demonstrated with results based on particle simulations, both for systems under 



temperature gradient [|10|,|TT| and shear ||12|| . 

Ordinarily the macroscopic study of subsonic hydrodynamical instabilities is based on 
the incompressibility assumption. However, as flrst pointed out by Zaitsev and Shliomis 
p!3| , this assumption is basically inconsistent with the very foundations of the fluctuating 



hydrodynamics formalism since it imposes flctitious correlations between the velocity com- 
ponents of the fluid. On the other hand, the compressibility of the fluid affects mostly fast 
sound modes, whereas the dynamics of the system near an instability is governed by slow 
dissipative modes. We may thus expect that the behavior of a fluid evolving near a sub- 
sonic instability threshold is practically not affected by its compressibility. This intuitive 
argument has been used by many authors who have considered fluctuating incompressible 
hydrodynamic equations, or even directly the corresponding normal form amplitude equa- 
tions to which they added random noise terms Jl^. In these approaches, the characteristics 
of the noise terms cannot be related to equilibrium statistical properties of the fluid and 
thus remain arbitrary. A more satisfactory approach would be to start with the full com- 
pressible fluctuating hydrodynamic equations. Reducing these equations to a flnal normal 
form amplitude equation near the instability, would lead directly to the explicit form of the 
associated noise terms consistent with such requirements as the fluctuation-dissipation theo- 
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rem. Such a procedure, however, proves to be quite difficult mainly because of the boundary 
conditions. To our knowledge, the only attempt in this direction has been made by Schmitz 



and Cohen for the case of Benard instability |[TH| . Concentrating on the behavior of a small 



layer in the bulk, these authors have succeeded in deriving the linearized fluctuating equa- 
tions close to the convective instability. Whether this technique can be generalized to derive 
the corresponding normal form amplitude equation for the case of Benard instability is not 
clear at present time. 

Recently, we have considered the problem of hydrodynamic fluctuations in the case of 
a simple flow proposed some fifty years ago by Kolmogorov ||16|. Thanks to the periodic 
boundary conditions associated to this model, a detailed analysis of the linearized fluctuating 
hydrodynamic equations, from near equilibrium up to the vicinity of the ffist instability, 
could be carried out In particular, we have been able to show that in the long time limit 
the flow behaves basically as if the fluid were incompressible, regardless of the value of the 
Reynolds number. The situation was different for the short time behavior. We established 
that the incompressibility assumption leads here to a wrong form of the static correlation 
functions, in agreement with Zaitsev and Shliomis prediction , except near the instability 



threshold, where our results strongly suggest that the incompressibility assumption becomes 
again valid. On the other hand, the linearized fluctuating hydrodynamic equations are 
clearly not valid close to, or beyond the instability threshold. Although extensive numerical 
simulations have basically confirmed our predictions, a satisfactory answer to this important 
problem requires a full nonlinear analysis of the fluctuating Kolmogorov flow. The present 
article is devoted to this problem. 

In the next section, the Kolmogorov flow is briefly reviewed. A nonlinear analysis is 
carried out for an incompressible fluid and the explicit form of the stream function and the 
associated velocity field is derived above but close to the instability. Section III is devoted 
to the analysis of a compressible fluid. After setting up a perturbation scheme, we show 
that the solution of the problem is basically the same as the one derived in section II for 
the incompressible fluid, at least close to the instability threshold. We then concentrate on 



the statistical properties of the flow and show that, close to the instability threshold, the 
dynamics of the system is governed by a set of two nonlinear coupled Langevin equations. 
Here again, the equivalence with the incompressible case is established. Concluding remarks 
and perspectives are summarized in section IV. 



II. INCOMPRESSIBLE KOLMOGOROV FLOW 

Consider an isothermal flow in a rectangular box Lj. x Ly oriented along the main axes, 
that is {0 < a; < L^,, < ?/ < Ly}. Periodic boundary conditions are assumed in both direc- 
tions and the flow is maintained through an external force field of the form 

Fea:i = -^0 sin (2 TT n y/Lj^) la; , (l) 

where 1^ is the unit vector in the x direction. This model represents the so-called Kolmogorov 
flow and it belongs to the wider class of two-dimensional negative eddy viscosity flows . 



It is entirely characterized through the strength of the force field Fq, the parameter n, which 
controls the wave number of the forcing, and the aspect ratio a^, defined as 

Oir — Lx/ Ly . (2) 

In the following, we will mainly concentrate on the case n = 1. 
The fluctuating hydrodynamic equations for this model read: 

V.(pv), (3) 

= -p(v V)v - Vp - V-<7 + Fe,t, (4) 

where p is the mass density, p the hydrostatic pressure and cr the two dzmensiona/ fluctuating 
stress tensor: 
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S is a random tensor whose elements {Sij} are Gaussian white noises with zero mean and 
covariances given by: 



< ^,,(r, t) S,,,{r', t') > = 2kBTo S{t - t') S{r - r') [r/(<5f;5f/ + ^f/^jj) + (C - vKP^J 



(6) 

For simphcity, we shall assume that the shear and bulk viscosity coefficients, t] and (, are 
state independent, i.e. they are constant. 

Let us first concentrate on the deterministic behavior. It can be easily checked that at 
the stationary state, the pressure and the density are uniform in space (p^t = Po, Pst = Po), 
whereas the velocity proffie is given by: 

Vst = uo sin (2 ny/Ly) 1^, 

For small enough Fq, this stationary flow is stable. As we increase Fq, however, the flow 
eventually becomes unstable giving rise to rotating convective patterns. Other instabilities 
of increasing complexity may occur for larger values of Fq, culminating in a chaotic like 
behavior similar to what is observed in turbulent flows 0-2O|. In this paper we shall limit 



ourselves to the analysis of the system near its first instability. 

We still have to supply the momentum conservation equation (Q) with an equation of 
state relating the pressure to the density (recall that the system is isothermal). In this 
section, we shall simply assume that the flow is incompressible, i.e. 

^ du dv ^ 

dx dy 

where u and v represent the x and y components of the velocity, respectively, i.e. v = 
ulx + vly. Relation (|^) implies a uniform density po throughout the system for all time, if 
initially so, as well as the existence of a scalar function ip{x,y), known as stream function, 
defined by the relations: 

u = ^, v = (9) 
dy ' dx 



+ R-^ (V^ ^p) + 8 TT^ R~' cos (2 n y) (10) 
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Scaling lengths by Ly, velocity by uq and time by Ly/uo, the dimensionless equation for the 
stream function reads: 

(9(V2^) _ d^diV^ijj) dipd{V'^i)) 
dt dy dx dx dy 

where R is the Reynolds number: 

R 

The stationary solution of ( p!OD reads: 



PO UQ Ly 

V 



(11) 



Z 71 



(12) 



Setting = ipst + Sip, and linearizing ([10|) around V'st, one gets 

^ ^ = - sm(27ry)^^— ^ - 47r2 sm(27r|/)^ + /^-^ ( 5^^) . (13) 
Owing to periodic boundary conditions, 6ilj{x,y,t) can be expanded in Fourier series: 

oo 

Sij{x,y,t)= exp{-2iTikyy) exp{-27iik^x/ar)SiJk,,ky(t) , 

/■I 1 

Si^k^,ky{t) = dy exp {2 n i ky y) — dx exp {2n i k^x/ar) Stp^x, y, t) . 

(14) 



Equation (|T^ can be then transformed to : 



dt 



where we have set 



-An'R-' (ki + ky^iP,. 



+ 7rk^ dipk^^ky + i - Sipk^^ky-i 

SlPk^^ky + l + SlPk^,ky-l 



+ 27r^^"^^ 



kl + ky 



(15) 



kx/Ojj, 



(16) 



In its general form, the analysis of this equation proves to be quite difficult pT|. On 
the other hand, if ipst is stable then, in the long time limit, the evolution of the system 
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will be mainly governed by long wavelength modes. Accordingly, we start our analysis by 
considering only the modes ky = 0,±1, i.e. we assume that Stplk^ , kyt) ^ for \ky\ > 2 
2^ . Defining the vector <5i/'fc^ = {S4'k^,o , ^i^k^,i ? ? the equation (|T5|) can be written 



in the following matricial form : 



dt 



(17) 



with 



A 



~47r^R-'kl 
7,K{l-kl)/{l + kl) 

y-^l^[i-ki)/[i + ki) 



\ 



-ATi^R-\l + kl) 

-Ati^R-\1 + kl 

We first note that the matrix A is diagonal for /cj, = so that the solution of eq. (|1 
simply reduces to: 



(18) 



5^0, i(i) ~ 5^0, -i(i) ~ exp(-47r^i?"'t). 



(19) 



Furthermore, by definition of the stream function, eq. (^, ipo^oit) = 0, Wt. We thus concen- 
trate on the case k^ ^ 0, looking for a similarity transformation T- A-T^^ which diagonalizes 
the matrix A. After some algebra, one finds 

/ 



(Ai - \^)/iik^ 



1 -1 



(A2-A3)/7r^, 1 



-1 



V 

where {Aj} are the eigenvalues of A: 







(20) 



Ai = -2t:^R~\1 + 2kl) + ^2kl{l - kl) / {I + kl) + A-n^ R-\ 
X, = -27r'R~'{l + 2~kl) - ^,j2kl{l-kl)/{l + kl) + A7,^R-\ 
A3 = -ATi'^R-^il + kl). 



(21) 



The equation ([T7|) then becomes: 

dt 



\iS(l)^it), I = 1,2,3 



(22) 



where 



d(f) = T ■ Si^ (23) 

It follows from pl| ) that A2 and A3 are always negative, whereas there exists a critical 
value of the Reynolds number 

Rcih) = 2^271 l±jL ; <~kl < 1 (24) 

'i - k 



for which Ai vanishes, thus indicating the limit of stability of the corresponding mode p3 
Clearly is an increasing function of {k^l, so that the first modes to become unstable 
correspond to \kx\ = 1, provided the aspect ratio > 1. As I, Rc —>■ 00, indicating 

that no instability can develop for perturbations of the same spatial periodicity as the applied 
force [Q. In the following, we shall therefore concentrate mainly on the case > 1. 



For = 2, relation ( p^ ) predicts a critical Reynolds number of Rc ~ 12.8255. Analytical 
calculations can still be handled when the modes ky = ±2 are taken into account as well, 
and lead to 

'1/2 



Rf\k^) = R,{k^) 



i.^x + 3) 



< kl < I (25) 



2{kl + Af{kl-l)\ 

For a.r = 2, one finds a critical Reynolds number of R^^'^ ~ 12.8738, so that the discrepancy 
remains below 0.4%. Numerical evaluation of Rc performed with a total amount of 103 modes 
shows no further significant discrepancy. We thus conclude that one can rely reasonably well 
on a "3-modes approximation theory" (that is Sipkx,ky{t) ~ for \ky\ > 2). It remains to 
check whether this approximation leads to the correct velocity field beyond the instability. 
To this end we need to work out the explicit form of the stream function. 

The calculations are tedious and quite lengthy, so that here we only report the basic 
steps. We start with the full nonlinear evolution equation for Sip = — tpst'- 

— sm[2 7ry) — An sm(z7ryj 



dt dx dx 

Jdtfj (9(V25^) d6ip d{ 

dy dx dx dy 
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+ fl- V^^Hn,) - ^ + ^ . (26) 



As for the linear case, we take the Fourier transform of this equation, hmiting ourselves 
to the first three modes ky = 0,±1. Applying then the transformation T to the resulting 
equation, cfr. eq. (|23|), one obtains: 

= A.50.(t) + , t = 1,2,3 , (27) 

where the $j's are nonlinear polynomial functions of 50i,502 and 503 and their complex- 
conjugates. Close to the bifurcation point {R ^ Rc , = I), the mode 6(f)i exhibits a 
critical slowing down since Ai(fc^ = 1) ^ 0. On this slow time scale, i.e. t ~ 0{X{^), the 
fast modes 6(f)2 and 503 can then be considered as stationary, their time dependence arising 
mainly through 501 (t). Setting d5(j)2/dt ~ 9503 /9t ~ 0, one can express the fast modes 
502 and 503 in terms of the slow mode, 50i, and its complex-conjugate, 50i*. If now one 
inserts the so-obtained expressions of the fast modes into the evolution equation of the slow 
mode, one obtains a closed nonlinear equation for the latter {adiabatic elimination [^,^). 
In practice, however, such a calculation is possible only close to the the bifurcation point, 
where the amplitude of 50i is supposed to approach zero as i? — > Rc. In fact, there exist 
other types of transitions, such as the one arising in the Vanderpol equation, where the 
amplitude of the solution above the instability does not vanish as one approaches to the 



critical point p^. Detailed analysis shows that this is not the case here (i.e. |50i| 
as i? — > -Rc), so that we can limit ourselves to lowest orders in |50i|, obtaining finally the 
so-called normal form or amplitude equation for the slow mode: 

-7|50i(t)r50i(t) (l + 0{\5Mt)\')) , (28) 

where 

^-Mk.^l)^'4 - I) O (W". - 11% (29) 

and 7 is a positive constant whose expression, to dominant order in \R/Rc — 1|, is given by 

/- 3 (g^ + I7at + 16a'r - 32) (a,^ + 1)^ 
^" a3 (a2 - 1)3/2 + 2)3 (a2 + 4)2 " W 
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Above the bifurcation point R > Rc {X > 0), the amphtude equation (pS]) admits two 
stable stationary solutions, corresponding to the rotation sense of the stream lines in the 
fluid: 



IX 
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= ±\l'- exp{i9o) 



(31) 



where 6o is a constant whose value depends on the initial conditions. The fact that the 
stationary solution still depends on the initial conditions simply reflects the Galilean in- 
variance in the x direction which results from the periodic boundary conditions imposed on 
the system. Using relation (^), one can compute the explicit form of the fast modes for 
= 0, ±1, ±2. Applying the inverse transform T~^{k^), cfr. eq. (pOf), to the so-obtained 
vector 6(j)^{kx) = {S(f)f, 6(f)f, and taking its inverse Fourier transform, one obtains the 
explicit expression of the stream function in real space. Up to order 0{R/Rc — 1), one gets: 



1 



^sti^: y) = -7^cos{2ny) ± 



271 
dr R( 



2 7r(a2 + 2) 
sin {2TTx/aj. — 6*0) sin(27r|/) 



+ 



Rl 



2ix{al + 2y 



C0s{2t{X / ttr — Oq) 



■ COS {Anx/ar — 20q 



cos(27ry) 



(32) 



where we have set = Using relations 

obtained straightforwardly: 



the velocity profiles can now be 



Usti.^, y) = sin(27r?/) t 
R^ 



An 



(a2 + 2) 



(a2 + 2y 



\S(f)i\ sin {2Txx/ar — 6q) cos(27r?/) 



cos {Anx/ar — 29q 



(a2 + 4)- 



sin(27r?/) 



(33) 



vfti^, y) = ± 



Rr 



[al + 2) 



sin {2nx/ar — Oq) + 



A-n 

CLr Rr 



COS {2nx/ar — Oq) sin(27r?/) 



2 R^ 



sin {Anx/ar — 20o) cos(27ry) 



(34) 



(a2 + 2)2 (a2 + 4)2 

A density plot of the stream function ( P^ is represented in Figure (1) for R = 15, ar = 2 
and ^0 = where, for the sake of clarity, a vector plot of the velocity field is also included. 
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We note that the flow has an ABC-hke topology with closed streamlines (eddies), open 
ones and separatrices between them. 

We recall that the above results rest on the "three modes" approximation theory. To 
check the validity of this basic assumption, we have solved numerically the incompressible 
nonlinear hydrodynamic equations for = 2, using standard techniques [^. Figure (2) 
compares contour plots of the stream function obtained numerically with its corresponding 
theoretical counterpart, eq. (^), for R = 15. Given the relatively large distance from 
the critical point {R/Rc — 1 ~ 17%), the agreement is much better than expected, the 
discrepancy remaining below 5%. Surprisingly, the agreement does not improve as we con- 
sider smaller values of the Reynolds number. This is shown in Figure (3), where both the 
numerical and the theoretical horizontal proflles of the stream function with a flxed value 
of the vertical coordinate, y = 3/4, are depicted for the Reynolds number R = 13. The 
discrepancy now exceeds 10%. 

To understand the origin of this unexpected behavior, we note that the value of the 
critical Reynolds number that we have used to evaluate the stream function (eq. |52|) is based 
on the three modes approximation theory (cfr. eq. ^). As shown before, the accuracy of 
the latter value of Rc is about 0.4%, which is flne as far as the distance from the critical 
point {R/Rc — 1) remains much larger than 0.4%. Now, for R = 13, the distance from the 
critical point is about 1% which is of the same order as the accuracy of Rc and explains the 
relatively important discrepancy we have observed in Figure (3). 

To overcome this difficulty, one has to compute a more accurate value of the critical 



Reynolds number, based for instance on the 5 modes approximation theory (cfr. eq. [25| ). As 
well known [Q, this correction concerns only the value of Rc, and in no way compromises 
the validity of the amplitude equation (^) and its corresponding solution eq. (p2D . This 
is illustrated in Figure (3), where excellent agreement with the numerical result is demon- 
strated, whenever we use R'^^ as the critical Reynolds number. For smaller value of R, one 
can as well compute numerically the value of Rc with desired precision and used it as an 
input to the amplitude equation (|28|) . 
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So far, we have limited ourselves to the analysis of the deterministic equations only, 
i.e. we have discarded the noise terms. In principle, there is no difficulty in taking into 
account the noise contributions as well, except that the amplitudes of the field variables 
(501,502,503) are now directly related to the amplitude B of the noise, which is typically 
a small parameter. For example, the fast variables (5025 503) ~ 0(i3^/^), whereas the slow 
variable 50i ~ 0(i3^/^) (a detailed discussion of this problem is given in fSlI]). Keeping this 
restriction in mind, one can repeat all the above calculations in the presence of noise terms. 
To the dominant order in |50i|, one finds 

^^^'^^^ = A50i(t) - 7|50i(t)p50i(t) + m, 



dt 

^^"^'^^^ A50Ut) -^\SMt)\'5m + at). (35) 



dt 

The functions ^{t) and its complex-conjugate ^*(t) are Gaussian white noises with zero 
means and correlations given by: 

<mm> =0 

<mCit')> =B5{t-t') (36) 



with 

_ 4 ke Tq al 
MulR 

M being the total mass of the system: 



1 + 0(|i?/i?,-l|) 



(37) 



M = arPoLl . (38) 

The results derived in this section were based explicitly on the incompressibility assump- 
tion. However, as discussed in the Introduction, this assumption is inconsistent with the 
very foundation of the fiuctuating hydrodynamic formalism. On the other hand, we have 
presented in numerical evidence that in the vicinity of the bifurcation point the system 
behaves basically as an incompressible fiuid. We therefore expect that the Langevin equa- 
tion ( ^51) should remain valid for R close enough to Re- We shall clarify this main issue in 
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the next section. 



III. FLUCTUATIONS IN THE COMPRESSIBLE FLOW 

Let us now consider the compressible hydrodynamic equations - ^ for which we need 
to specify an equation of state. Since the system is isothermal, we simply set 

P = cIp, (39) 

where is the isothermal speed of sound. As in the previous section, we start with the 
linearized hydrodynamic equations around the reference state {pQ.Wst}, where v^j is given 
by eq.(0). Setting 

p = Pq + 5p, 

V = + 5v , (40) 

and scaling lengths by Ly, time by Ly/cs, Sp by po and 5v by the speed of sound Cg, the 
dimensionless linear fluctuating equations in the Fourier space read (recall that kx = kx/ar): 

dSpk^,ky{t) ^2vri (kxSuk,,ky + kySvk^^ky) + e Rti K {Spk,,ky+i - Spk,,ky-i) , (41) 

^^^y ^ + Svk,,ky~i) + ne Rkx{Suk^,ky+i - ^Uk,,ky-i) 

- 4:7T^e{kl + kl)5uk,,ky - 4:7r'^ ae~kx{kxSuk,,ky + ky6vk,,ky) 
+ 2TTikx6pk,,ky + Fk^ , ky (t) , (42) 

^^^''g^' =7TeR~kx (fafc., ky+i - Svk^,ky-i) - 4:71'^ e{kl + kl)5vk^,ky 

- A7r'^aeky{kx5uk,,ky + ky^^k^,ky) + 2 7r i ky 5pk,,ky + Gk^,ky{t) , (43) 



where R is the Reynolds number, defined in eq.([TT 

e = ^ 
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(44) 



and 

a = Ch. (45) 

The functions F^^^k^ and Gk^^ky are Fourier components of the noise terms ; their covariances 
follow directly from eqs. (||, ||): 

< Fu^^ky it) Gfc. , (f) > =8n^eAa K ky 5^+^,0 ^it - t') , 
<G,,,.,(t)G,,,,,(t') > = 8n'eA[kl + {a + 1) kl]S^;^,^,6{t - t') , 



where k = {kx, ky) and 



(46) 



Mc,2' ^^^^ 



M being the total mass of the system (cfr. eq. (|38|) ). 

For the sake of clarity, we first focus on the deterministic behavior, i.e. we discard for 
the moment the noise contributions from the evolution equations (0- 43). Furthermore, we 



shall limit ourselves to the 3-modes approximation theory, i.e. we shall neglect the modes 
with \ky\ > 2, for the very same reasons that we have discussed for the incompressible case. 



With these assumptions, eqs. (^ - reduce to a system of nine coupled equations. It can 
then be checked that the change of variables 

Spi{t) = 5pk^,i{t) ±5pk^,^i{t) 

Svlit) = 6v,^,,it)±6v,^,^,it) (48) 

leads to a "partial diagonalization" of the evolution equations, i.e. the equations for the 
variables {5pfc^,o, (^Pfc^? '^'^fc^,0; ^'^k^y ^"^fc^} decouple from the rest. Furthermore, their asso- 
ciated eigenvalues prove to remain strictly negative, regardless of the value of the Reynolds 
number R, so that they are not determinant for the onset of convective instability. We 
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therefore focus on the remaining four variables {^Sp^^Su'^^, Svj^^^, Svk^^Y Defining the vector 
^^kAt) = {Spt, Sut^, Svg, 6vk,fi} one readily finds: 

d 



dt 



(49) 



where the matrix C is given by: 

I 



27iikx 

2Tci —An'^eakx 




27ri 

—An'^eakx 
Ati^I + a + kl) 



irekxR 





-27ieR 
-27reKR 



-A-n^kl 



(50) 



The analysis can be simplified somewhat by noticing that the parameter e must remain 



small if one wishes to remain within the limit of validity of the hydrodynamic regime [32 



Furthermore, as already mentioned in the introduction, in this article we limit ourselves to 
strictly sub-sonic flows, so that we shall restrict the analysis to a parameter domain where 



e -C 1 ; eR = uq/cs -C 1. 



(51) 



Accordingly, we evaluate the eigenvalues of the matrix C perturbatively: 



~X{k,) = ~X^^\k,) + £A«(fc,) + 



After some algebra, one finds, up to order 0(£:^): 



Ai(fc^ 



~ 2 



-2 7r^(l + 2k'J + 7rV47r2 + 2 R^ k, (1 - kl)/{l + k^) 



~ 2 



-2vr^(l + 2ki) - rrdAn^ + 2 R^ k^ (1 - kl)/{l + k^) 



\s{k,) = 2711^1 + kl- 2'K^{a + 1)^(1 + kl), 
~X^(k^) = -2Txi Jl + kl - 2 7r2(a + 1) e{l + kl) . 



(52) 



(53) 
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The eigenvalues Ai and A2 correspond to dissipative (viscous) modes, while A3 and A4 are 
related to the propagation of (damped) sound waves. It can then be easily checked that the 
real parts of A2, A3 and A4 are always negative, whereas there exists a critical value of the 
Reynolds number 

lip 

RciK) = 2V2 7r , ^ ; < fc^ < 1 , (54) 
^1 - ~kl 



for which Ai vanishes, thus indicating the limit of stability of the corresponding mode. 

Remarkably, the above expression of the critical Reynolds number is identical to the 
one obtained in the incompressible case (cfr. eq.(p^)). In fact, detailed analysis shows that 
the relation ( p^ is exact, i.e. it is independent of e, at least within the framework of the 
3-modes approximation theory. On the other hand, if the modes ky = ±2 are taken into 
account as well, tedious calculations lead to 

-1/2 



R^^ikx) — Rc{kx 



_l_ ^x i^x + 3) 



+ O {{uo/csy) ; < < 1 (55) 



2ikl + Ayikl-l)_ 

which is again equivalent to the corresponding result obtained for the incompressible case, 
eq. (|25|), the correction being of the order of O(e^). In particular, the first mode to become 
unstable corresponds to \kx\ = 1, provided a.^ > 1. 

We note that the matrix C is singular for kx = 0, i.e. one of its eigenvalues vanishes. 
A close inspection shows that this zero eigenvalue corresponds to the mode 5fo,o which is 
identically zero because of linear momentum conservation. Accordingly, in what follows we 
shall concentrate on the case k^^ ^ 0, looking for a similarity transformation S ■ C ■ 
which diagonalizes the matrix C. For consistency, here again we perform the calculations 
perturbatively, i.e. we expand S in powers of e: 

Sik,) = So(M + eS^ik,) + ... (56) 



Note that this method constitutes an alternative to the time scale perturbation theory ||3^ 
that was generalized by Schmitz and Cohen ||T3 in order to study the Benard instability in 



a compressible fluid. 
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Since the explicit form of the eigenvalues are known up to O(e^), we only need to eval- 
uate S (and its inverse S^^) up to the same order. Despite this simplification, the general 
expression of S is quite awkward and will not be presented here. The rest of the calculations 
are quite straightforward, but remain tedious and lengthy, so we only give a brief sketch of 
the basic steps (see the discussions below eq. (^) ). 

We start by taking the Fourier transform of the fluctuating hydrodynamic equations 
(§ - H). Using the change of variables (^0|) and (|48|) , we next derive the nonlinear fluctuating 
equations for ^h^^ . We then apply the transformation S to the latter, obtaining a set of four 
nonlinear equations for the variables 502, 503, ^04} = S4>{t) = S ■ 6hk^. Close to the 
bifurcation point {R ^ Rc , = 1), the mode 50i exhibits a critical slowing down, since 
by construction Ai ~ 0. We can therefore proceed to an adiabatic eliminations of the "fast" 
modes |502, 503, 504|, limiting ourselves to dominant orders in |50i| (see the paragraph 
preceding eq. (|35|) ) . The final result is a set of two coupled Langevin equations for the slow 
mode 501 and its complex conjugate 50^ : 

~i{t) 



dt 
dt 



A50i(t) - 7|50i(t)|'50i(t) + m 
~X64>m - 7|50i(t)p50tW +C{t) 



(57) 



with 



A = Ai(A;,. = l) = A ^ 



and 



7 = 7 — 

Mo 



1 + 0{ul/c'^ 



2 + l 

4:71 e ■ 



aj{aj + 2) 



7 

eR 



(58) 



(59) 



where A and 7 are given by eqs. (pQ]) and (pOD, respectively. The functions ^(t) and its 
complex-conjugate ^*(t) are Gaussian white noises with zero means and correlations given 
by: 



<i{t)m> =0, 

<lm\^)> =B6{t-t') 



(60) 



with 

3 



1 + 0{ul/cl)\ ^ AealA, (61) 

where B and A are given by eqs. ( p^ ) and (^7D, respectively. 

Although the form of the Langevin equations ( |57D is the same as the one obtained for the 
incompressible case, eqs. (pSf), they are nevertheless not equivalent since their coefficients are 
clearly different, even to dominant order in e. The main reason for this apparent discrepancy 
is related to the fact that, for the incompressible case, the analysis has been carried out by 
scaling the velocities by mq, whereas for the compressible case we have used a different scaling, 
i.e. we have scaled the velocities by the velocity of sound Cg. If now we switch back to the 
former scaling, i.e. we perform the change of variables t ^ t Cg/uo, {u,v} uo/cs{u,v}, 
then eqs. (^) lead to 

6Mt) = -SMt) [l + 0{ul/cl)] (62) 

Remarkably, this result shows that, to dominant order in e, the evolution of fluctuating 
compressible and incompressible hydrodynamic equations are governed by the very same 
slow mode, at least for values of Reynolds number close to its critical value. 

Let us first consider the macroscopic behavior. Using eqs. (p5D, ( |5D| ) and (p2|), one can 
go backward step by step and derive as well the evolution equations of the hydrodynamical 
velocities near the instability threshold. It can then be easily checked that, to dominant 
order in e, the compressible stationary velocity profiles are given by their incompressible 
expressions, eqs. (|33| and 0). To check this important result, we have solved numerically the 
full nonlinear compressible hydrodynamic equations and compared the result with analytical 
expressions obtained for the incompressible case. A typical result is shown in Figure (4), 
where Ust{x,y = 1/4) as a function of Vst{x,y = 1/4) is depicted for R = 15, e = 10^^ and 

= 2. Given the relatively large values of the Reynolds number {R/Rc — 1 ~ 17%) and e, 
the agreement is very good, the discrepancy remaining below 5%. 

We now concentrate on the behavior of fiuctuations, as described by the Langevin equa- 



tions (1571) . The associated Fokker-Planck equation reads: 
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dt 



d 



d 



-(A 501 - 750?(50t)P + 



dm) 

At the stationary state, one finds: 



1 «0f ) 



P + 



B dP 

2dm) 

B dP 

2 dm 



(63) 



P. 



N exp 



lixm? -lm\' 



(64) 



with 



AT = i ^Ti B/^ exp (aV 7i3) erfc (-\/^^B 



(65) 



where erfc(. . .) stands for the complementary error function. Thanks to this result, one 
readily gets: 



< |50i|2 >= 4 f A + B /AM 



7 

Away from the bifurcation point (A << 0) the quartic term in 
the distribution is Gaussian and 



(66) 

is negligible so that 
(67) 



2 7r2(i?2_/?2) (^2 + 1) 

The fluctuations thus behave as ~ 0{A^^'^). Recall that the parameter A is inversely 
proportional to the system's total number of particles so that ^ << 1 (cfr. eq. (^7D ). As 
one approaches the bifurcation point, the Gaussian character of the distribution is gradually 
lost. Right at the bifurcation point, A = 0, one has 



< m\'^ >A=o = 2ear 



1/2 



(68) 



which shows that the fluctuations now behave as ~ (9(^^/^). The enhancement of 

fluctuations and the change of the probability law at the bifurcation point are a direct 
manifestation of spatial symmetry breaking associated with the emergence of convective 
patterns. 
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On the other hand, the fast modes |502, (^03, S^a\ prove to remain Gaussian, regardless 
of the value of the Reynolds number. Detailed analysis shows that their contribution to 
nonequilibrium statistical properties of the fluid remain of the order of 0{uq/c1). In other 
words, the fluctuation spectrum of hydrodynamic variables are mainly determined by the 
statistical properties of 6 (pi. For instance, the static velocity auto-correlation function is 
found to obey: 



< 5vk ■ 5v_k > -2A= 4y^^ < > 



1 + O Uuo/cs) 



(69) 



a2 (a2 + 2)2 

where the second term on the left hand side is the equilibrium contribution and < > 
is given by eq. (|66D . 

It is instructive to study the Gaussian limit, R << Rc, where the linearized Langevin 



equations, eqs. (|41| - |^), remain valid. As has been shown in they lead to the following 



expression for the static velocity auto-correlation function: 

< . iv_. >,-2A= (70) 

Now, inserting into eq. ( |50D the Gaussian form of < >, as given by eq. (§3), leads 

precisely to the very same result. We thus conclude that our general expression, eq. (pUD, 
remains valid in the Gaussian regime, R << Rc, despite the fact that it has been derived in 
the close vicinity of the bifurcation point R ~ Rc- 

To check the validity of our theoretical results, we have simulated the nonlinear fluctu- 
ating hydrodynamic equations - ^) for different values of R, setting 
and A = 10~V256 ^ 3.9 x 10"^ The estimated statistical errors remain below 5% for 
R < 10, but grows rapidly as we consider higher values of R, reaching about 13% for 
R ^ Rc- Above the bifurcation point, R > Rc, the stationary distribution has two maxima. 



located at 6(f)i = ±yX/^, which correspond (up to a phase factor) to the deterministic 
stationary solutions of the amplitude equation (pTj). Because of the presence of noise terms, 
the system visits these states in a rather random fashion, resulting into a huge dispersion 
of data. This is specially true for R close to Rc, which is precisely the situation where our 
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theoretical predictions are expected to be applicable. Under this circumstance, obtaining 
reliable statistics requires prohibitively large computing times, so that we have been forced 
to limit the numerical simulations to values of Reynolds numbers R < Re- 

The results are presented in Figure (5), together with both the complete and the lin- 
earized solutions, eqs. and (|TDD, respectively. The linear theory (Gaussian limit) shows 
quantitative agreement for values of R/Rc up to about 86%, but significant discrepancies 
start to show up as i? — > Rc where the theory leads to diverging correlation functions (cfr. 
eq. ([70|) ). This is not the case for the complete solution, eq. (|69|) , which exhibits perfect 
quantitative agreement for R/Rc up to 95%. A relatively small discrepancy of about 8% is 
however observed for higher values of R. Although this discrepancy remains within the limit 
of the estimated statistical errors, its systematic aspect requires nevertheless some clarifica- 
tions. In this respect, it is important to recall that the results derived in this section were 
valid up to 0(mq/c^). Now, by definition uq/cs = Re (cfr. eq. (pi])), and since we have 
set e = 10^^, Rc6 ^ 0.13 at the bifurcation point. This relatively large value of RcE might 
well be at the origin of the observed discrepancy. To check the validity of this argument, 
it is tempting to perform the simulations all over again for a smaller value of e. However, 
since the relaxation time of hydrodynamical modes grows as e~^, reaching the same degree 
of statistical accuracy as for the previous cases requires much more longer running times. 
For this reason we decided to perform only one more simulation right at the critical point, 
R = Rc, setting e = 10^^. The theoretical prediction for the nonequilibrium part of the 
velocity correlation function is 2.31 x 10~^. The simulation leads to 2.24 x 10^® with an 
estimated statistical errors of about 15%. The discrepancy is now about 3%, much better 
than for the case e = 10^^. 
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IV. CONCLUDING REMARKS 



Recently, we have studied the statistical properties of the linearized Kolmogorov flow, 
from near equilibrium up to the vicinity of the first instability leading to the formation of 
vortices [0. In particular, we have established that the incompressibility assumption leads 
to a wrong form of the static correlation functions, except near the instability threshold 
where numerical results suggest that the incompressibility assumption should remain valid. 
The clarification of this important issue requires a nonlinear analysis of the fluctuating 
Kolmogorov flow. This is precisely the main purpose of the present article. 

We have first considered the case of an incompressible fluid. After identifying the slow 
modes, governing the evolution of the system in the vicinity of the instability threshold, we 
have performed an adiabatic elimination of the fast modes to obtain a set of two nonlinear 
Langevin equations for the slow modes. We have then succeeded to derive the explicit 
form of the stationary stream function, as well as the corresponding velocity profiles, in real 
space. Numerical studies of the nonlinear hydrodynamical equations allowed us to confirm 
our theoretical predictions. 

We have next considered the case of the compressible Kolmogorov flow. The analysis can 
be simplified somewhat by noticing that the evolution of a compressible fluid is generally 
characterized by two different time scales: a slow one, related to the dissipative viscous 
modes, and a fast one, expressing the propagation of (damped) sound modes. The ratio of 
these time scales, denoted by e (cfr. eq. ^H), can be considered as a small parameter, since 
otherwise the very validity of the hydrodynamics can no more be guaranteed [^. We thus 
have at our disposal a natural small parameter which can be used to set up a perturbative 
technique. As already mentioned, this method constitutes an alternative to the time scale 
perturbation theory that was generalized by Schmitz and Cohen in order to study the Benard 
instability in a compressible fluid ||T5| , |33| . 

Using this perturbation technique, we have first shown that the macroscopic behavior 
of the fluid is not affected - up to C(mq/c^) - by the compressibility, in agreement with 
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the intuitive arguments that we have presented in the Introduction. We then succeeded to 
estabhsh that, close to the instabihty threshold, the stochastic dynamics of the system is 
governed by two coupled non linear Langevin equations in the Fourier space. The solution 
of these equations can be cast into the exponential of a Landau Ginzburg functional which, 
to dominant order in £, proves to be identical to the one obtained for the case of the incom- 
pressible fluid. The theoretical predictions have been confirmed by numerical simulations of 
the nonlinear fluctuating hydrodynamic equations. 



ACKNOWLEDGMENTS 

We are very grateful to professors E. G. D. Cohen, G. Nicolis, J. W. Turner and C. Van 
den Broeck for helpful comments. This work is supported by the Belgian Federal Office for 
Scientific, Technical and Cultural Affairs within the framework of the "Poles d'attractions 
interuniversitaires" program, and by an European Commission DC 12 Grant PSS*1045. 



24 



REFERENCES 

[1] I. Bena, M. Malek Mansour, and F. Baras, Phys. Rev. E 59, 5503 (1999). 

[2] L. D. Landau and E. M. Lifshitz, Fluid Mechanics, Pergamon, Oxford (1984). 

[3] T. R. Kirkpatrick, E. G. D. Cohen and J. R. Dorfman, Phys. Rev. Lett. 42, 862 (1979) ; 
ibid. 44, 472 (1980) ; ibid. Phys. Rev. A 26, 995 (1982). 

[4] J. Dufty, in Spectral Line Shapes, P. Wende ed., De Kruger, Berhn (1981) ; J. Lutsko 
and J. Dufty, Phys.Rev. A 32, 1229 (1985). 

[5] L Procaccia, D. Ronis, I. Oppenheim, Phys. Rev. Lett. 42, 287 (1979) ; D. Ronis, I. 
Procaccia and L Oppenheim, Phys. Rev. A 19, 1324 (1979) ; A. -M. Tremblay, M. Arai 
and E. Siggia, Phys. Rev. A 23, 1451 (1981) ; D. Ronis, L Procaccia and L Oppenheim, 
Phys. Rev. A 26, 1812 (1982). 

[6] G. Van der Zwan, D. Bedeaux and P. Mazur, Physica A 107, 491 (1981) ; R. Schmitz, 
and E.G.D. Cohen, J.Stat.Phys. 39, 285 (1985) ; ibid. 40, 431 (1985). 

[7] J. Machta, L Procaccia and L Oppenheim, Phys. Rev. Lett. 42, 1368 (1979) ; J. Machta, 
I. Oppenheim and I. Procaccia, Phys. Rev. A 22, 2809 (1980). 

[8] R. Schmitz, Phys. Reports 171, 1 (1988). 

[9] B.M. Law, and J.V. Sengers, J.Stat.Phys. 57, 531 (1989); 

B.M. Law, P.N. Segre, R.W. Gammon, and J.V. Sengers, Phys.Rev. A 41, 816 (1990). 

[10] M. Malek Mansour, A.L. Garcia, G.Lie, and E. Clementi, Phys. Rev. Lett. 58, 874 
(1987) ; M. Mareschal, M. Malek Mansour, G. Sonino, and E. Kestemont, Phys.Rev. A 
45, 7180 (1992). 

[11] A. Suarez, J. P. Boon and P. Grosfils, Phys. Rev. E 54, 1208 (1996). 

[12] A.L. Garcia, M. Malek Mansour, G. Lie, M. Mareschal, and E. Clementi, Phys.Rev. A 
36, 4348 (1987). 

25 



[13] V.M. Zaitsev, and M.I. Shliomis, Sov.Phys. JETP 32, 866 (1971). 

[14] R. Graham, Phys. Rev. A 10, 1762 (1974) ; R. Graham and H. Pleiner, Phys. Fluids, 
18, 130 (1975); J. Swift and P. C. Hohenberg, Phys. Rev. A 15, 319 (1977). 

[15] R. Schmitz and E. G. D. Cohen, J. Stat. Phys. 39, 285 (1985); ibid. 40, 431 (1985). 

[16] A.M. Obukhov, Russ. Math. Survey 38,113 (1983). 

[17] S. Gama, M. Vergassola, and U. Frisch, J. Fluid Mech. 260, 95 (1994). 

[18] E. N. Lorenz, J. Fluid Mech. 55, 545 (1972). 

[19] Z. S. She, Phys. Lett. A 124, 161 (1987); ibidem. Ph. D. Thesis, Univ. Paris VII (1987); 
G. I. Sivashinsky, Physica D 17, 243 (1985). 

[20] R. Benzi, and S. Sued, J. Stat. Phys. 56, 69 (1989); U. Frisch, B. Legras, and B. Villone, 
Physica D 94, 36 (1996). 

[21] L. Meshalkin, and Ya. G. Sinai, J. Appl. Mah. Mech. (P.M.M.) 25, 1140 (1961). 

[22] J.S.A. Green, J. Fluid. Mech. 62, 273 (1974). 

[23] This result is not new and has been obtained by several authors, like for instance A. A. 
Nepomniaschichii, J. Appl. Math. Mech. (P.M.M.) 40, 836 (1976), Meshalkin and Sinai 
12111, and Lorenz m. 



[24] C. Machioro, Commun. Math. Phys. 105, 99 (1986). 

[25] H. Haken, Synergetics, An Introduction, 3rd ed.. Springer Ser. Synergetics, Vol. 1, 
Springer, Berlin, 1983. 

[26] G. Nicolis, Introduction to Nonlinear Science, Cambridge University Press, Cambridge 
(1995). 

[27] F. Baras, M. Malek Mansour, and C. Van Den Broeck, J. Stat. Phys. 28, 577 (1982). 



26 



[28] G. M. Zaslavsky, R. Z. Sagdeev, D. A. Usikov, and A. A. Chernikov, Weak chaos and 
quasi-regular patterns, Cambridge University Press, Cambridge (1991), chapter 9. 

[29] C.-Y. Chow, An Introduction to Computational Fluid Mechanics, Seminole PubL, Boul- 
der CO (1983). 

[30] This result has already been obtained by M. A. Brutyan and P.L. Krapivsky, Eur. J. 
Mech. Fluids B, 11, 587 (1992) for a compressible fluid in the limit of inflnite aspect 
ratio. 

[31] C. Van den Broeck, M. Malek Mansour, and F. Baras, J. Stat. Phys. 28, 557 (1982). 

[32] See for instance W. E. Alley and B. J. Alder, Phys. Rev. A 27, 3158(1983) ; W. E. 
Alley, B. J. Alder and S. Yip, Phys. Rev. A 27, 3174 (1983). 

[33] U. Geigenmiiller, U. M. Titulaer and B. U. Felderhof, Physica 119 A, 41 (1983). 



27 



FIGURES 

FIG. 1. Density plot of the stream function, eq. (|3^), for R = 15, a,. = 2 and = 0. For the 
sake of clarity, a vector plot of the velocity field is also included. 

FIG. 2. Stationary state contour plot of the stream function for R = 15, = 2 and = 0. 
The full and dashed lines correspond to theoretical prediction (eq. ^2p and numerical results, 
respectively. The discrepancy remains below 5%. 

FIG. 3. Horizontal profile of the stationary state stream function, with y = 3/4, as a function 
of the vertical coordinate x for = 13, = 2 and 9q = 0. The full and dashed lines represent 
theoretical predictions obtained by using an estimation of the critical Reynolds number based 
on 5-modes (eq. |2^) and 3-modes (eq. ^) approximation theories, respectively. The diamonds 
correspond to numerical results. 

FIG. 4. Vertical versus horizontal components of the stationary state velocity field with y = 3/4. 
The full line corresponds to theoretical predictions, as given by eqs. ( |33|) and (|3^), whereas the 
dashed line is obtained by solving numerically the compressible nonlinear hydrodynamic equations. 
The parameters are R = 15, = 2, 6*0 = and e = 10^^. The discrepancy is about 5%. 

FIG. 5. Fourier transform of the nonequilibrium part of the static velocity auto-correlation 
function, normalized by the corresponding equilibrium part, as a function of R/Rc- The solid and 
dashed curves represent the complete and the linearized solutions, eqs. ( p9| ) and (pO]), respectively, 
whereas the diamonds correspond to numerical results obtained by the simulation of nonlinear 
compressible fluctuating hydrodynamic equations. The parameters are = 2, e = 10^2 and 
A = 10~^/256. The estimated statistical error is about 13% for the last data point. 
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